Zadatak 1

U timu predmeta na MSTeams platformi dostupna je baza slikanih šaka koje pokazuju „papir“, „kamen“ i „makaze“. Projektovati inovativni sistem za prepoznavanje pokazanih znakova zasnovan na testiranju hipoteza.

a) Detaljno opisati algoritam za obradu slike i odabir obeležja koji prethodi samoj klasifikaciji. Algoritam treba da bude što robusniji (na različite osvetljaje, položaje šaka, načine pokazivanja znakova itd).

b) Izvršiti podelu na trening i test skup. Rezultate klasifikacije test skupa prikazati u obliku matrice konfuzije.

c) Odabrati dva znaka i dva obeležja takva da su odabrani znakovi što separabilniji u tom prostoru. Prikazati histogram obeležja za oba slova i prokomentarisati njihov oblik.

d) Za slova i obeležja pod c) projektovati parametarski klasifikator po izboru i iscrtati klasifikacionu liniju.

a)

Preprocesiranje slike

Za prepoznavanje znakova rukom (papir, kamen, makaze) na slikama sa zelenom pozadinom, potrebno je izvršiti niz koraka preprocesiranja kako bi se od originalne slike dobila binarna maska koja jasno izdvaja šaku.

Funkcija izoluj_saku kao ulaz prima putanju do slike, a vraća izdvojenu binarnu sliku šake.

Slika se prvo učitava u RGB formatu, a zatim se konvertuje u HSV (hue, saturation, value) format. Ovaj format je pogodniji za segmentaciju boja jer omogućava da se boje razdvoje nezavisno od osvetljenja i senki.

Pozadina slike je zelena, pa se definiše opseg nijansi zelene boje u HSV prostoru - donja_granica i gornja_granica. U tom opsegu izdvajaju se pikseli pozadine.

Određivanjem piksela koji nisu u definisanom opsegu zelene boje dobijamo masku šake. Saka ce biti bele boje, dok ce pozadina biti crna. Dobijena slika predstavlja binarnu (crno-belu sliku).

Radi bolje obrade, uklanja se po 5 piksela sa gornje, donje i leve ivice slike. Ovim se eliminišu potencijalno loše oblasti na granicama. Ne uzimamo i desnu ivicu zato sto se cesto tu nalazi šaka.

Da bi se popunile praznine i otklonile nepravilnosti u konturi, koristi se zatvaranje slike (closing):

• dilatacija popunjava rupe u segmentisanim delovima
• erozija nakon toga vraća konturu u granice sličnim početnim (istanjuje ono što je dilatacija podebljala)

Na dobijenu sliku primenjuje se median filter sa maskom $7 \times 7$. Na ovaj način se uklanja „salt & pepper“ šum koji se može javiti u procesu segmentacije.

Kao poslednji korak, slika se kropuje tako da ostane samo deo gde se nalazi šaka. Leva granica se pomera nadesno dok ne naiđe na kolonu sa belim pikselima. Sličan posao obavljam za ostale ivice.

Na ovaj način se dobija kropova binarna slika šake, koju ću koristiti za klasifikaciju.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import train_test_split
from scipy.stats import mode
from skimage import io, color, morphology, filters
import os

def izoluj_saku(putanja_slike):
    
    img = io.imread(putanja_slike)
    hsv = color.rgb2hsv(img)
    donja_zelena = np.array([85/360, 0.2, 0.2]) 
    gornja_zelena = np.array([170/360, 1.0, 1.0])
    
    maska_pozadine = (
        (hsv[..., 0] >= donja_zelena[0]) & (hsv[..., 0] <= gornja_zelena[0]) &
        (hsv[..., 1] >= donja_zelena[1]) & (hsv[..., 1] <= gornja_zelena[1]) &
        (hsv[..., 2] >= donja_zelena[2]) & (hsv[..., 2] <= gornja_zelena[2])
    )
    
    maska_sake = ~maska_pozadine
    
    binarna = np.where(maska_sake, 255, 0).astype(np.uint8)
    binarna = binarna[5:binarna.shape[0]-5, 5:]
    
    dilatacija = morphology.dilation(binarna, morphology.disk(2))
    erozija = morphology.erosion(dilatacija, morphology.disk(2))

    zatvaranje = erozija.astype(np.uint8) * 255
    
    filtrirana = filters.median(zatvaranje, np.ones((7, 7), dtype=np.uint8)) * 255
    
    M,N = np.shape(filtrirana)

    left = 0
    while(left < N-1 and np.sum(filtrirana[:, left] == 255) == 0):
        left += 1
    
    right = N-1
    while(right > 0 and np.sum(filtrirana[:, right] == 255) == 0):
        right -= 1
    
    up = 0
    while(up < M-1 and np.sum(filtrirana[up, :] == 255) == 0):
        up += 1
    
    down = M-1
    while(down > 0 and np.sum(filtrirana[down, :] == 255) == 0):
        down -= 1

    kropovana = filtrirana[up:down, left:right]

    return kropovana

rock_processed = [] 
scissors_processed = []
paper_processed = []

path = r'E:\SEDMI SEMESTAR\PO'
folder_path = os.path.join(path, 'RPS')
folder_paper = os.path.join(folder_path,'paper')
folder_rock = os.path.join(folder_path,'rock')
folder_scissors = os.path.join(folder_path,'scissors')

for filename in os.listdir(folder_scissors):
    image_path = os.path.join(folder_scissors, filename)
    scissors_processed.append(izoluj_saku(image_path))

for filename in os.listdir(folder_rock):
    image_path = os.path.join(folder_rock, filename)
    rock_processed.append(izoluj_saku(image_path))
    
for filename in os.listdir(folder_paper):
    image_path = os.path.join(folder_paper, filename)
    paper_processed.append(izoluj_saku(image_path))

Primer

image_path = os.path.join(folder_scissors, os.listdir(folder_scissors)[2])
image = izoluj_saku(image_path)

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow(io.imread(image_path))
ax[0].set_title("Neprocesirana slika (makaze)")
ax[0].axis("off")

ax[1].imshow(image, cmap='gray')
ax[1].set_title("Procesirana slika (makaze)")
ax[1].axis("off")

plt.tight_layout()
plt.show()


image_path = os.path.join(folder_rock, os.listdir(folder_rock)[2])
image = izoluj_saku(image_path)

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow(io.imread(image_path))
ax[0].set_title("Neprocesirana slika (kamen)")
ax[0].axis("off")

ax[1].imshow(image, cmap='gray')
ax[1].set_title("Procesirana slika (kamen)")
ax[1].axis("off")

plt.tight_layout()
plt.show()


image_path = os.path.join(folder_paper, os.listdir(folder_paper)[0])
image = izoluj_saku(image_path)

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow(io.imread(image_path))
ax[0].set_title("Neprocesirana slika (papir)")
ax[0].axis("off")

ax[1].imshow(image, cmap='gray')
ax[1].set_title("Procesirana slika (papir)")
ax[1].axis("off")

plt.tight_layout()
plt.show()

Atributi

Za atribute sam uzela udeo belih piksela u levom delu slike, tj. gledala sam od šestine do trećine slike i u tom delu sračunala udeo belih piksela. Očekujem da makaze imaju najmanji udeo, dok bi kamen trebalo da ima najveći udeo, a papir će biti negde u sredini.

Za drugi atribut je sračunat centar šake, očekujem da kamen ima najveće vrednosti - jer je relativno centrirana šaka, dok preostale dve bi trebalo da imaju nešto manje vrednosti, jer je kod njih šaka ispružena (izduženija, nije u centru).

def udeo_belih(image):
    sirina = image.shape[1]
    
    start = int(np.floor(sirina / 6))
    stop = int(np.floor(sirina / 3))
    
    levi_deo = image[:, start:stop]
    
    beli_pikseli = np.sum(levi_deo == 255)
    ukupno_pikseli = levi_deo.size
    
    return beli_pikseli / ukupno_pikseli

def nadji_centar(image):
    rows, cols = np.where(image == 255)
    
    r = int(np.mean(rows))
    c = int(np.mean(cols))
    
    return r, c

def izdvoj_atribute(images):
    atributi = np.zeros((2, len(images)))
    
    for i, image in enumerate(images):
        atributi[0, i] = udeo_belih(image)
        r, c = nadji_centar(image)
        atributi[1, i] = r/c
    
    return atributi

b)

features_paper = izdvoj_atribute(paper_processed)
features_rock = izdvoj_atribute(rock_processed)
features_scissors = izdvoj_atribute(scissors_processed)

s_dim = features_scissors.shape[1]
r_dim = features_rock.shape[1]

X = np.concatenate((features_scissors, features_rock, features_paper), axis=1)
Y = np.ones((X.shape[1]))
Y[s_dim:s_dim+r_dim] = 2
Y[s_dim+r_dim:] = 3

X_train, X_test, Y_train, Y_test = train_test_split(X.T, Y, test_size=0.2, stratify=Y, random_state=25)
X_train = X_train.T
X_test = X_test.T
X1 = X_train[:,Y_train==1]
X2 = X_train[:,Y_train==2]
X3 = X_train[:,Y_train==3]

plt.figure(figsize=(10,8))
plt.plot(X1[0,:], X1[1,:], 'x')
plt.plot(X2[0,:], X2[1,:], 'x')
plt.plot(X3[0,:], X3[1,:], 'x')
plt.legend(['makaze', 'kamen', 'papir'])
plt.title('Prikaz atributa')
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()

M1 = np.mean(X1, axis=1, keepdims=True)
M2 = np.mean(X2, axis=1, keepdims=True)
M3 = np.mean(X3, axis=1, keepdims=True)

S1 = np.cov(X1)
S2 = np.cov(X2)
S3 = np.cov(X3)

def gaussian_pdf(X, M, S):
    detS = np.linalg.det(S)
    invS = np.linalg.inv(S)
    n = S.shape[1]
    f = 1/(2*np.pi*detS**0.5)*np.exp(-0.5*(X-M).T@invS@(X-M))
    return f[0,0]


def plot_cm(Y, Y_pred, skup="trening"):
    cm = confusion_matrix(Y, Y_pred)

    plt.figure(figsize=(6, 5))
    sns.heatmap(cm, annot=True, fmt="d")

    plt.xlabel("Predikcija")
    plt.ylabel("Istina")
    plt.title(f"Konfuziona matrica za {skup} skup")
    plt.show()
    
    
def predict(X, M_list, S_list):
    Y_pred = np.zeros(X.shape[1])
    for i in range(X.shape[1]):
        odb = X[:, i].reshape(-1, 1)
        fs = []
        for M, S in zip(M_list, S_list):
            fs.append(gaussian_pdf(odb, M, S))
        Y_pred[i] = np.argmax(fs) + 1
    return Y_pred

Y_train_pred = predict(X_train, [M1, M2, M3], [S1, S2, S3])
plot_cm(Y_train, Y_train_pred)

Y_test_pred = predict(X_test, [M1, M2, M3], [S1, S2, S3])
plot_cm(Y_test, Y_test_pred, 'test')

Dobija se visok procenat tačnosti i na 3 klase koristeći samo 2 atributa

c)

Očigledno je da su kamen i makaze najrazdvojeniji tako da sam se odlučila da klasifikujem ta dva znaka. Pošto su linarno separabilni, najlakše je projektovati linearni klasifikator.

X = np.concatenate((features_scissors, features_rock), axis=1)
Y = np.zeros((X.shape[1]))
Y[features_s.shape[1]:] = 1

plt.figure()
sns.countplot(x=Y)
plt.title('Udeo podataka')
plt.show()

X_train, X_test, Y_train, Y_test = train_test_split(X.T, Y, test_size=0.2, stratify=Y, random_state=42)
X_train = X_train.T
X_test = X_test.T
X1 = X_train[:,Y_train==0]
X2 = X_train[:,Y_train==1]

plt.figure(figsize=(10,8))
plt.plot(X1[0,:], X1[1,:], 'x')
plt.plot(X2[0,:], X2[1,:], 'x')
plt.legend(['makaze', 'kamen'])
plt.title('Prikaz atributa')
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()

plt.figure()
plt.hist2d(X1[0,:], X1[1,:], bins=20, alpha=1, cmap='Blues')
plt.title('Histogram klase scissors')
plt.xlabel('atribut 1')
plt.ylabel('atribut 2')
plt.show()

plt.figure()
plt.hist2d(X2[0,:], X2[1,:], bins=20, alpha=1, cmap='Reds')
plt.title('Histogram klase rock')
plt.xlabel('atribut 1')
plt.ylabel('atribut 2')
plt.show()

Uočavamo da su centroidi značajno udaljeni

plt.hist(X1[0,:])
plt.title('Histogram klase scissors: Atribut 1')
plt.show()

plt.hist(X1[1,:])
plt.title('Histogram klase scissors: Atribut 2')
plt.show()

plt.hist(X2[0,:])
plt.title('Histogram klase rock: Atribut 1')
plt.show()

plt.hist(X2[1,:])
plt.title('Histogram klase rock: Atribut 2')
plt.show()

def druga_num_metoda(X1, X2):
    N1 = X1.shape[1]
    N2 = X2.shape[1]
    
    M1_est = np.mean(X1, axis=1).reshape(-1,1)
    M2_est = np.mean(X2, axis=1).reshape(-1,1)
    S1_est = np.cov(X1)
    S2_est = np.cov(X2)    
    
    s = np.arange(0,1,10**-3)
    v0_opt_s = []
    Neps_s = []
    
    for i in range(len(s)):
        V = np.linalg.inv(s[i]*S1_est + (1-s[i])*S2_est)@(M2_est-M1_est)
        Y1 = V.T@X1
        Y2 = V.T@X2
        Y = np.sort(np.concatenate((Y1, Y2), axis=1)) 
        v0 = np.zeros(Y.shape[1]-1) 
        Neps = np.zeros(Y.shape[1]-1)
        for j in range(Y.shape[1]-1):
            v0[j] = -(Y[0,j]+Y[0,j+1])/2
            Neps[j] += np.sum(Y1>-v0[j])
            Neps[j] += np.sum(Y2<-v0[j])  
        Neps_s.append(np.min(Neps))
        v0_opt_s.append(v0[np.argmin(Neps)])
    Neps_opt = min(Neps_s)
    v0_opt = v0_opt_s[np.argmin(Neps_s)]
    s_opt = s[np.argmin(Neps_s)]
    return (s_opt, v0_opt, Neps_opt, M1_est, M2_est, S1_est, S2_est)

s_opt, v0_opt, Neps_opt, M1_est, M2_est, S1_est, S2_est = druga_num_metoda(X1, X2)
V = np.linalg.inv(s_opt*S1_est + (1-s_opt)*S2_est)@(M2_est-M1_est)
v0 = v0_opt

x1 = np.arange(0.3,1,0.01)
x2 = -(v0+V[0]*x1)/V[1]

plt.figure(figsize=(10, 8))
plt.plot(X1[0,:],X1[1,:], 'x')
plt.plot(X2[0,:],X2[1,:], 'x')
plt.plot(x1,x2,'m')
plt.title('Klasifikacija klasa Makaze i Kamen')
plt.legend(['Makaze', 'Kamen'])
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()

Y_pred = (V.T@X_train+v0)
Y_pred = Y_pred.reshape(-1,)
Y_pred[Y_pred>=0] = 1
Y_pred[Y_pred<0] = 0

plot_cm(Y_train, Y_pred)

Y_pred = (V.T@X_test+v0)
Y_pred = Y_pred.reshape(-1,)
Y_pred[Y_pred>=0] = 1
Y_pred[Y_pred<0] = 0

plot_cm(Y_test, Y_pred, 'test')

Dobijamo tačost od 100% na oba skupa - linearno separabilni, pa linearni klasifikator odrađuje dobar posao.

Zadatak 2

Generisati po $𝑁 = 500$ odbiraka iz dveju dvodimenzionih bimodalnih klasa: $$\Omega_1 \sim P_{11} \cdot N(M_{11}, \Sigma_{11}) + P_{12} \cdot N(M_{12}, \Sigma_{12})$$ $$\Omega_2 \sim P_{21} \cdot N(M_{21}, \Sigma_{21}) + P_{22} \cdot N(M_{22}, \Sigma_{22})$$

Parametre klasa samostalno izabrati.

a) Na dijagramu prikazati odbirke.

b) Iscrtati kako teorijski izgledaju funkcije gustine verovatnoće za raspodele klasa i uporediti ih sa histogramom generisanih odbiraka.

c) Projektovati Bajesov klasifikator minimalne greške i na dijagramu, zajedno sa odbircima, skicirati klasifikacionu liniju. Uporediti grešku klasifikacije konkretnih odbiraka sa teorijskom greškom klasifikacije prve i druge vrste za datu postavku.

d) Projektovati klasifikator minimalne cene tako da se više penalizuje pogrešna klasifikacija odbiraka iz prve klase.

e) Ponoviti prethodnu tačku za Neuman-Pearson-ov klasifikator. Obrazložiti izbor $\epsilon_2 = \epsilon_0$.

f) Za klase oblika generisanih u prethodnim tačkama, projektovati Wald-ov sekvencijalni test pa skicirati zavisnost broja potrebnih odbiraka od usvojene verovatnoće grešaka prvog, odnosno drugog tipa.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from sklearn.metrics import confusion_matrix

import seaborn as sns
import scipy
from skimage import io
import os
from sklearn.metrics import confusion_matrix

a)

Prvo podešavamo parametre za dati problem

np.random.seed(0)
N = 500 
P1, P2 = 0.5, 0.5

P11, P12 = 0.3, 0.7
M11 = np.array([-0.7, 0.6]).reshape(-1, 1)
S11 = np.array([[0.7, 0.2], [0.2, 0.7]])
M12 = np.array([4, 4]).reshape(-1, 1)
S12 = np.array([[2, -0.3], [-0.3, 2]])

P21, P22 = 0.4, 0.6
M21 = np.array([-3, 1]).reshape(-1, 1)
S21 = np.array([[1, 0.7], [0.7, 1]])
M22 = np.array([-5, 4]).reshape(-1, 1)
S22 = np.array([[1, 0.1], [0.1, 1]])

Zatim vršimo nasumično uzorkovanje iz odgovarajućih raspodela - ovo predstavlja generisanje naših odbiraka.

X11 = np.random.multivariate_normal(M11.reshape(-1),S11,N)
X12 = np.random.multivariate_normal(M12.reshape(-1),S12,N)
X21 = np.random.multivariate_normal(M21.reshape(-1),S21,N)
X22 = np.random.multivariate_normal(M22.reshape(-1),S22,N)

pp = np.random.rand(1,N)
p = np.concatenate((pp,pp),axis=0).T
X1 = (p<P11)*X11 + (p>=P11)*X12

pp = np.random.rand(1,N)
p = np.concatenate((pp,pp),axis=0).T
X2 = (p<P21)*X21 + (p>=P21)*X22

Prikaz odbiraka

plt.figure()
plt.plot(X1[:,0],X1[:,1],'rx')
plt.plot(X2[:,0],X2[:,1],'yx')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['Klasa1', 'Klasa2'])
plt.show()

Postoji preklapanje između klasa - zbog ovoga klasifikatori neće imati potpunu tačnost, ali je ovakav primer podoban za ilustraciju karakteristika različitih klasifikatora

b)

Prikaz fgv

x = np.arange(-10,10.1,0.1)      #levu granicu inkluzivno, desnu eksluzivno
y = np.arange(-10,10.1,0.1)

f11 = np.zeros((len(x),len(y)))
f12 = np.zeros((len(x),len(y)))
f21 = np.zeros((len(x),len(y)))
f22 = np.zeros((len(x),len(y)))
f1 = np.zeros((len(x),len(y)))
f2 = np.zeros((len(x),len(y)))

def gaussian_pdf(X, M, S):
    detS = np.linalg.det(S)
    invS = np.linalg.inv(S)
    f = 1/(2*np.pi*detS**0.5)*np.exp(-0.5*(X-M).T@invS@(X-M))
    return f[0,0]

for i in range(0,len(x)):
    for j in range(0,len(y)):
        X = np.array([x[i],y[j]]).reshape(-1,1)
        f11[i,j] = gaussian_pdf(X,M11,S11)
        f12[i,j] = gaussian_pdf(X,M12,S12)
        f21[i,j] = gaussian_pdf(X,M21,S21)
        f22[i,j] = gaussian_pdf(X,M22,S22)
        f1[i,j] = P11*f11[i,j] + P12*f12[i,j]
        f2[i,j] = P21*f21[i,j] + P22*f22[i,j]

X, Y = np.meshgrid(x,y)    

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(X,Y,f1.T,cmap='viridis')
ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_title('Funkcija gustine verovatnoće f1 (bimodalna)')
plt.show()

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(X,Y,f2.T,cmap='viridis')
ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_title('Funkcija gustine verovatnoće f2 (bimodalna)')
ax.view_init(elev=30, azim=45)
plt.show()

class1 = np.vstack((X11, X12))
class2 = np.vstack((X21, X22))

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.hist2d(class1[:, 0], class1[:, 1], bins=30, cmap='Blues')
plt.colorbar(label='Broj tačaka')
plt.contour(X,Y,f1.T)
plt.title("Klasa 1 - 2D histogram")
plt.xlabel("Prva dimenzija")
plt.ylabel("Druga dimenzija")

plt.subplot(1, 2, 2)
plt.hist2d(class2[:, 0], class2[:, 1], bins=30, cmap='Reds')
plt.colorbar(label='Broj tačaka')
plt.contour(X,Y,f2.T)
plt.title("Klasa 2 - 2D histogram")
plt.xlabel("Prva dimenzija")
plt.ylabel("Druga dimenzija")

plt.tight_layout()
plt.show()

Jasno se vidi preklapanje izohipsi fgv sa hitogramom generisanih odbiraka. Ovo potvrđuje tačnost koda.

c) Bajesov test minimalne verovatnoće greške

Pretpostavimo dve klase $\omega_1$ i $\omega_2$. Za svaki odbirak (tj. za svaki oblik) potrebno nam je znanje o:

• apriornim verovatnoćama pojava klasa $P(\omega_1)=P_1$ i $P(\omega_2)=P_2$
• gustinama verovatnoće uslovljenim klasom $f(x\mid \omega_1)=f_1(x)$ i $f(x\mid \omega_2)=f_2(x)$.

U praksi $f_i(x)$ često modelujemo Gaussovom raspodelom, a vektor srednje vrednosti i kovariacionu matricu određujemo na osnovu apriornog znanja ili na osnovu odbiraka.

Klasifikator ima za zadatak da na osnovu datih elemenata sračuna aposteriorne verovatnoće za obe klase, uporedi sa granicom i donese odluku.

$\frac{f_1(x)}{f_2(x)}$ poredimo sa pragom $\frac{P_2}{P_1}$. Ukoliko je veći od praga, odbirak pripada klasi 1, u suprotnom pripada klasi 2.

Klasifikaciona linija se često zapisuje i kao: $$h(x) = -ln \left( \frac{f_1(x)}{f_2(x)} \right)$$ tada se prag svodi na $$T = ln \left( \frac{P_1}{P_2} \right)$$ Zbog promene znaka ovde odbirak pripada klasi 1 ako je diskriminaciona funkcija manja od praga, a pripada klasi 2 ukoliko je veća.

U ovom delu uporedićemo pristup gde znamo fgv i koristi parametre koje smo ranije definisali i pristup gde aproksimiramo parametre i koristimo (unimodalnu) gausovsku fgv. Ovo radim kako bih pokazala da je adekvatno nekada modelovati različite fgv gausovskim.

h = np.zeros((len(x),len(y)))
T = np.log(P1/P2)

K1 = np.empty((0,2))
K2 = np.empty((0,2))

for i in range(0,len(x)):
    for j in range(0,len(y)):
        X = np.array([x[i],y[j]]).reshape(1,2)
        h[i,j] = -np.log(f1[i,j]/f2[i,j])
        if h[i,j] < T:
            K1 = np.append(K1,X,axis=0)
        else:
            K2 = np.append(K2,X,axis=0)
            
plt.figure()
plt.plot(X1[:,0],X1[:,1],'rx')
plt.plot(X2[:,0],X2[:,1],'yx')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['Klasa1', 'Klasa2'])
plt.scatter(K1[:, 0], K1[:, 1], color='lightcoral', alpha=0.3, s=10, edgecolor='none')  
plt.scatter(K2[:, 0], K2[:, 1], color='palegreen', alpha=0.3, s=10, edgecolor='none')
plt.contour(x, y, h.T, levels=[T], colors='black')
plt.title('Klasifikaciona kriva sa pravom fgv')
plt.show()

data = np.concatenate((X1,X2),axis=0)

y_true = np.zeros(X1.shape[0]+X2.shape[0]) 
y_true[X1.shape[0]:] = 1
y_predict = np.zeros(np.shape(y_true))

for i in range(0,len(y_true)):
    sample = data[i,:].reshape((-1, 1))
    h = -P11 * np.log(gaussian_pdf(sample,M11,S11)) - P12 * np.log(gaussian_pdf(sample,M12,S12)) + P21 * np.log(gaussian_pdf(sample,M21,S21)) + P22 * np.log(gaussian_pdf(sample,M22,S22))
    y_predict[i] = h>T

cm = confusion_matrix(y_true,y_predict)
plt.figure()
sns.heatmap(cm, annot=True, fmt='g')
plt.title('Konfuziona matrica za pravu fgv')
plt.show()

Sračunaćemo greške dobijene ovom metodom

eps1b = cm[0,1]/sum(cm[0,:]) 
eps2b = cm[1,0]/sum(cm[1,:])
print(f'Greška prve vrste iznosi: {eps1b:.3f}')
print(f'Greška prve vrste iznosi: {eps2b:.3f}')
e1 = 0
e2 = 0

for i in range(0,len(x)):
    for j in range(0,len(y)):
        if -np.log(f1[i,j]/f2[i,j]) < 0:
            e2 += f2[i,j]*0.1*0.1
        else:
            e1 += f1[i,j]*0.1*0.1
            
print('Greška prve vrste teorijski iznosi: ', e1)
print('Greška druge vrste teorijski iznosi: ', e2)

Greška prve vrste iznosi: 0.036
Greška prve vrste iznosi: 0.010
Greška prve vrste teorijski iznosi:  0.020083858465993084
Greška druge vrste teorijski iznosi:  0.01513947433866279

M1 = np.mean(X1, axis=0).reshape(-1,1)
M2 = np.mean(X2, axis=0).reshape(-1,1)
S1 = np.cov(X1.T)
S2 = np.cov(X2.T)

h = np.zeros((len(x),len(y)))
T = np.log(P1/P2)

K1 = np.empty((0,2))
K2 = np.empty((0,2))

for i in range(0,len(x)):
    for j in range(0,len(y)):
        X = np.array([x[i],y[j]]).reshape(1,2)
        f1[i,j] = gaussian_pdf(X.T,M1,S1)
        f2[i,j] = gaussian_pdf(X.T,M2,S2)
        h[i,j] = -np.log(f1[i,j]/f2[i,j])
        if h[i,j] < T:
            K1 = np.append(K1,X,axis=0)
        else:
            K2 = np.append(K2,X,axis=0)
            
plt.figure()
plt.plot(X1[:,0],X1[:,1],'rx')
plt.plot(X2[:,0],X2[:,1],'yx')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['Klasa1', 'Klasa2'])
plt.scatter(K1[:, 0], K1[:, 1], color='lightcoral', alpha=0.3, s=10, edgecolor='none')#scatter, a ne plot da se ne bi povezale sve tacne i bio bag     
plt.scatter(K2[:, 0], K2[:, 1], color='palegreen', alpha=0.3, s=10, edgecolor='none')
plt.contour(x, y, h.T, levels=[T], colors='black')
plt.title('Klasifikaciona kriva sa estimiranom fgv')
plt.show()

data = np.concatenate((X1,X2),axis=0)
M1 = np.mean(X1, axis=0).reshape(-1,1)
M2 = np.mean(X2, axis=0).reshape(-1,1)
S1 = np.cov(X1.T)
S2 = np.cov(X2.T)

y_true = np.zeros(X1.shape[0]+X2.shape[0])
y_true[X1.shape[0]:] = 1
y_predict = np.zeros(np.shape(y_true))
for i in range(0,len(y_true)):
    sample = data[i,:].reshape((-1, 1))
    h = -np.log(gaussian_pdf(sample,M1,S1))+np.log(gaussian_pdf(sample,M2,S2))
    y_predict[i] = h>T

cm = confusion_matrix(y_true,y_predict)
plt.figure()
sns.heatmap(cm, annot=True, fmt='g')
plt.title('Konfuziona matrica za estimiranu fgv')
plt.show()

Uočavamo da klasifikaciona kriva koja koristi stvarnu fgv daje bolje rezultate, me]utim i klasifikaciona kriva koja koristi procenjenu fgv odra]uje dobar posao. Takođe uočavamo da klasifikaciona kriva sa procenjenom fgv predstavlja krivu drugog reda, dok je ona koja koristi pravu fgv kriva višeg reda.

Sada računamo greške dobijene na osnovu procenjene fgv

eps1b = cm[0,1]/sum(cm[0,:]) 
eps2b = cm[1,0]/sum(cm[1,:])
print(f'Greška prve vrste iznosi: {eps1b:.3f}')
print(f'Greška prve vrste iznosi: {eps2b:.3f}')
e1 = 0
e2 = 0

for i in range(0,len(x)):
    for j in range(0,len(y)):
        if -np.log(f1[i,j]/f2[i,j]) < 0:
            e2 += f2[i,j]*0.1*0.1
        else:
            e1 += f1[i,j]*0.1*0.1
            
print('Greška prve vrste teorijski iznosi: ', e1)
print('Greška druge vrste teorijski iznosi: ', e2)

Greška prve vrste iznosi: 0.050
Greška prve vrste iznosi: 0.012
Greška prve vrste teorijski iznosi:  0.04787945243437909
Greška druge vrste teorijski iznosi:  0.027052146602492902

Razlika u greškama se javlja zbog greške prilikom numeričke integracije, kao i zbog relativno malog broja odbiraka nad kojima procenjujemo eksperimentalnu grešku.

d) Bajesov test minimalne cene

Proces je analogan procesu projektovanja Bajesovog klasifikatora minimalne verovatnoće greške. Ovde dolazi samo do jednog dodatka, uvodimo cene, gde će cena $C_{ij}$ predstavljati cenu odluke da odbirak pripada klasi i, a on u stvari pripada klasi j. Sada umesto očekivane verovatnoće greške minimizujemo očekivani rizik, koji uračunava cene $$r = min(r_1(x), r_2(x))$$ $$r_1(x) = c_{11}q_1(x) + c_{12}q_2(x)$$ $$r_2(x) = c_{21}q_1(x) + c_{22}q_2(x)$$

Oblik klasifikacione funkcije je isti kao kod Bajesovog klasifikatora minimalne verovatnoće greške, jedino se menja izraz za prag, koji sada postaje:$$T = ln \left( \frac{C_{21} - C_{11}}{C_{12} - C_{22}} \cdot \frac{P_1}{P_2} \right)$$

Prema tekstu zadatka treba da važi $C_{21} > C_{12}$.

c11 = 0
c22 = 0
c12 = 1
c21 = 4
T = np.log(P1/P2*(c21-c11)/(c12-c22))

h = np.zeros((len(x),len(y)))

for i in range(0,len(x)):
    for j in range(0,len(y)):
        X = np.array([x[i],y[j]]).reshape(-1,1)
        f11[i,j] = gaussian_pdf(X,M11,S11)
        f12[i,j] = gaussian_pdf(X,M12,S12)
        f21[i,j] = gaussian_pdf(X,M21,S21)
        f22[i,j] = gaussian_pdf(X,M22,S22)
        f1[i,j] = P11*f11[i,j] + P12*f12[i,j]
        f2[i,j] = P21*f21[i,j] + P22*f22[i,j]

K1 = np.empty((0,2))
K2 = np.empty((0,2))

for i in range(0,len(x)):
    for j in range(0,len(y)):
        X = np.array([x[i],y[j]]).reshape(1,2)
        h[i,j] = -np.log(f1[i,j]/f2[i,j])
        if h[i,j] < T:
            K1 = np.append(K1,X,axis=0)
        else:
            K2 = np.append(K2,X,axis=0)
            
plt.figure()
plt.plot(X1[:,0],X1[:,1],'rx')
plt.plot(X2[:,0],X2[:,1],'yx')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['Klasa1', 'Klasa2'])
plt.scatter(K1[:, 0], K1[:, 1], color='lightcoral', alpha=0.3, s=10, edgecolor='none')
plt.scatter(K2[:, 0], K2[:, 1], color='palegreen', alpha=0.3, s=10, edgecolor='none')
plt.contour(x, y, h.T, levels=[T], colors='black')
plt.show()

data = np.concatenate((X1,X2),axis=0)

y_true = np.zeros(X1.shape[0]+X2.shape[0]) 
y_true[X1.shape[0]:] = 1
y_predict = np.zeros(np.shape(y_true))

for i in range(0,len(y_true)):
    sample = data[i,:].reshape((-1, 1))
    h = -P11 * np.log(gaussian_pdf(sample,M11,S11)) - P12 * np.log(gaussian_pdf(sample,M12,S12)) + P21 * np.log(gaussian_pdf(sample,M21,S21)) + P22 * np.log(gaussian_pdf(sample,M22,S22))
    y_predict[i] = h>T

cm = confusion_matrix(y_true,y_predict)
plt.figure()
sns.heatmap(cm, annot=True, fmt='g')
plt.title('Konfuziona matrica za pravu fgv')
plt.show()

eps1c = cm[0,1]/sum(cm[0,:]) 
eps2c = cm[1,0]/sum(cm[1,:])

print(f'Greška prve vrste iznosi: {eps1c:.3f}')
print(f'Greška druge vrste iznosi: {eps2c:.3f}')

Greška prve vrste iznosi: 0.018
Greška druge vrste iznosi: 0.016

Zaista uočavamo da je greška prve vrste sada manja, međutim zato je došlo do poraste greške druge vrste. Ponovićemo i za pravu fgv.

M1 = np.mean(X1, axis=0).reshape(-1,1)
M2 = np.mean(X2, axis=0).reshape(-1,1)
S1 = np.cov(X1.T)
S2 = np.cov(X2.T)


c11 = 0
c22 = 0
c12 = 1
c21 = 4
T = np.log(P1/P2*(c21-c11)/(c12-c22))

h = np.zeros((len(x),len(y)))

K1 = np.empty((0,2))
K2 = np.empty((0,2))

for i in range(0,len(x)):
    for j in range(0,len(y)):
        X = np.array([x[i],y[j]]).reshape(1,2)
        f1[i,j] = gaussian_pdf(X.T,M1,S1)
        f2[i,j] = gaussian_pdf(X.T,M2,S2)
        h[i,j] = -np.log(f1[i,j]/f2[i,j])
        if h[i,j] < T:
            K1 = np.append(K1,X,axis=0)
        else:
            K2 = np.append(K2,X,axis=0)
            
plt.figure()
plt.plot(X1[:,0],X1[:,1],'rx')
plt.plot(X2[:,0],X2[:,1],'yx')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['Klasa1', 'Klasa2'])
plt.scatter(K1[:, 0], K1[:, 1], color='lightcoral', alpha=0.3, s=10, edgecolor='none')#scatter, a ne plot da se ne bi povezale sve tacne i bio bag     
plt.scatter(K2[:, 0], K2[:, 1], color='palegreen', alpha=0.3, s=10, edgecolor='none')
plt.contour(x, y, h.T, levels=[T], colors='black')
plt.show()

data = np.concatenate((X1,X2),axis=0)
M1 = np.mean(X1, axis=0).reshape(-1,1)
M2 = np.mean(X2, axis=0).reshape(-1,1)
S1 = np.cov(X1.T)
S2 = np.cov(X2.T)

y_true = np.zeros(X1.shape[0]+X2.shape[0])
y_true[X1.shape[0]:] = 1
y_predict = np.zeros(np.shape(y_true))
for i in range(0,len(y_true)):
    sample = data[i,:].reshape((-1, 1))
    h = -np.log(gaussian_pdf(sample,M1,S1))+np.log(gaussian_pdf(sample,M2,S2))
    y_predict[i] = h>T

cm = confusion_matrix(y_true,y_predict)
plt.figure()
sns.heatmap(cm, annot=True, fmt='g')
plt.show()

eps1c = cm[0,1]/sum(cm[0,:]) 
eps2c = cm[1,0]/sum(cm[1,:])

print(f'Greška prve vrste iznosi: {eps1c:.3f}')
print(f'Greška druge vrste iznosi: {eps2c:.3f}')

Greška prve vrste iznosi: 0.004
Greška druge vrste iznosi: 0.074

e) Neyman-Pearson-ov test

Potrebno je odrediti $\epsilon_0$ i $\mu$. Prag je sada definisan kao: $$T = -ln~\mu$$

$\epsilon_0$ možemo inicijalizovati na zbir grešaka prvog i drugog tipa koje smo dobili za Bajesov klasifikator minimalne verovatnoće greške i nadati se da $\epsilon_1$ teži nuli, jer smo "sve greške" smestili u grešku drugog tipa, ali to ne mora uvek biti slučaj. Svakako ovo je dobra polazna tačka i dodatne promene $\epsilon_0$ mogu se vršiti nakon provere performansi klasifiktora.

Pošto nije jednostavno iz zadatog $\epsilon_0$ dobiti $\mu$, mi ćemo sračunati procenu greške drugog tipa ta različite vrednosti $\mu$, a potom odabrati $\mu$ za traženo $\epsilon_0$. Potom ćemo projektovati klasifikator.

Koristićemo estimiranu fgv

M1 = np.mean(X1, axis=0).reshape(-1,1)
M2 = np.mean(X2, axis=0).reshape(-1,1)
S1 = np.cov(X1.T)
S2 = np.cov(X2.T)

x = np.arange(-10,10.1,0.1)     
y = np.arange(-10,10.1,0.1)

f1 = np.zeros((len(x), len(y)))
f2 = np.zeros((len(x), len(y)))
h = np.zeros((len(x), len(y)))

for i in range(len(x)):
    for j in range(len(y)):
        X = np.array([x[i], y[j]]).reshape(2,1)
        f1[i,j] = gaussian_pdf(X,M1,S1)
        f2[i,j] = gaussian_pdf(X,M2,S2)
        h[i,j] = -np.log(f1[i,j]/f2[i,j])

mu = np.arange(0.01, 10, 0.01)
eps2 = np.zeros(len(mu))

for m in range(len(mu)):
    for i in range(len(x)-1):
        for j in range(len(y)-1):
            if h[i, j] < -np.log(mu[m]):
                eps2[m] = eps2[m] + (f2[i,j] + f2[i+1,j] + f2[i, j+1] + f2[i+1,j+1])/4*0.1*0.1                               

plt.figure()
plt.plot(mu,eps2)
plt.xlabel('Parametar $\mu$')
plt.ylabel('Greška druge vrste')
plt.show()    

eps0 = eps1b + eps2b
for i in range(len(eps2)):
    if(eps2[i]<=eps0):
        break
        
mu_opt = mu[i]
print('Optimalna vrednost ', mu_opt)

Optimalna vrednost  0.38

T = -np.log(mu_opt)

h = np.zeros((len(x),len(y)))

K1 = np.empty((0,2))
K2 = np.empty((0,2))

for i in range(0,len(x)):
    for j in range(0,len(y)):
        X = np.array([x[i],y[j]]).reshape(1,2)
        f1[i,j] = gaussian_pdf(X.T,M1,S1)
        f2[i,j] = gaussian_pdf(X.T,M2,S2)
        h[i,j] = -np.log(f1[i,j]/f2[i,j])
        if h[i,j] < T:
            K1 = np.append(K1,X,axis=0)
        else:
            K2 = np.append(K2,X,axis=0)
            
plt.figure()
plt.plot(X1[:,0],X1[:,1],'rx')
plt.plot(X2[:,0],X2[:,1],'yx')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['Klasa1', 'Klasa2'])
plt.scatter(K1[:, 0], K1[:, 1], color='lightcoral', alpha=0.3, s=10, edgecolor='none')#scatter, a ne plot da se ne bi povezale sve tacne i bio bag     
plt.scatter(K2[:, 0], K2[:, 1], color='palegreen', alpha=0.3, s=10, edgecolor='none')
plt.contour(x, y, h.T, levels=[T], colors='black')
plt.title('Neyman-Pearson-ov klasifikator')
plt.show()

data = np.concatenate((X1,X2),axis=0)
M1 = np.mean(X1, axis=0).reshape(-1,1)
M2 = np.mean(X2, axis=0).reshape(-1,1)
S1 = np.cov(X1.T)
S2 = np.cov(X2.T)

y_true = np.zeros(X1.shape[0]+X2.shape[0])
y_true[X1.shape[0]:] = 1
y_predict = np.zeros(np.shape(y_true))
for i in range(0,len(y_true)):
    sample = data[i,:].reshape((-1, 1))
    h = -np.log(gaussian_pdf(sample,M1,S1))+np.log(gaussian_pdf(sample,M2,S2))
    y_predict[i] = h>T

cm = confusion_matrix(y_true,y_predict)
plt.figure()
sns.heatmap(cm, annot=True, fmt='g')
plt.title('Konfuziona matrica za Neyman-Pearson-ov klasifikator')
plt.show()

I ovde se može uočiti da je greška prve vrste manja u odnosu na onu inicijalnu kod BTMVG.

eps1n = cm[0,1]/sum(cm[0,:]) 
eps2n = cm[1,0]/sum(cm[1,:])

print(f'Greška prve vrste iznosi: {eps1n:.3f}')
print(f'Greška druge vrste iznosi: {eps2n:.3f}')

Greška prve vrste iznosi: 0.012
Greška druge vrste iznosi: 0.058

f) Wald-ov test

Definišemo $S_m$ kao $$S_m = \sum_{i=1}^{m} h(X_i)$$, gde je $$h(x) = -ln \left( \frac{f_1(x)}{f_2(x)} \right)$$

Odluka se donosi na sledeći način $$Sm \left\{\begin{array}{lc}>b & X_i \in \omega_2 \\ < a & X_i \in \omega_1 \\ \in (a,~b) & \text {uzmi (m+1) odbirak} \end{array}\right.$$

$~~$ $$a = - ln \left( \frac{1-\epsilon_1}{\epsilon_2} \right)$$ $~$

$$b = - ln \left( \frac{\epsilon_1}{1-\epsilon_2} \right)$$

Ovaj test je pogodno koristiti kada nam odbirci stižu jedan po jedan (sekvencijalno) i kada treba donositi odluku u realnom vremenu.

Takođe za očekivanti broj potrebnih odbiraka prve, odnosno druge klase se pokazuje da važi: $$m_1 = \frac{(1-\epsilon_1)ln \left(\frac{\epsilon_2}{1-\epsilon_1} \right) - \epsilon_1 ln\left(\frac{\epsilon_1}{1-\epsilon_2} \right)}{\eta_1}$$ izraz za $m_2$ se dobnija analogno zamenom a i b u izraz $$m_2 = \frac{b(1-\epsilon_2) + a\epsilon_2}{\eta_2}$$

eps1w_arr = np.logspace(-7, -0.1, num=12)
eps2w_arr = np.logspace(-7, -0.1, num=12)
m1 = np.zeros((len(eps1w_arr), len(eps2w_arr)))
m2 = np.zeros((len(eps1w_arr), len(eps2w_arr)))
data = np.concatenate((X1,X2),axis=0).T

for i in range(len(eps1w_arr)):
    for j in range(len(eps2w_arr)):
        eps1w = eps1w_arr[i]
        eps2w = eps2w_arr[j]
        A = (1-eps1w) / eps2w
        B = eps1w / (1-eps2w)
        a = -np.log(A)
        b = -np.log(B)
        
        Sm = 0
        cnt = 0
        
        samples1 = []
        samples2 = []
        
        for n in range(data.shape[1]):
            f1 = gaussian_pdf(data[:,n].reshape(-1,1), M1, S1)
            f2 = gaussian_pdf(data[:,n].reshape(-1,1), M2, S2)
            h = -np.log(f1/f2)
            Sm += h
            cnt += 1
            
            if Sm < a:
                samples1.append(cnt)
                Sm = 0
                cnt = 0
            if Sm > b:
                samples2.append(cnt)
                Sm = 0
                cnt = 0
        
        if len(samples1) > 0:
            m1[i,j] = sum(samples1)/len(samples1)
        if len(samples2) > 0:
            m2[i,j] = sum(samples2)/len(samples2)  

plt.figure(figsize=(10,8))
for j in range(len(eps2w_arr)):
    plt.plot(eps1w_arr, m1[:, j])
plt.legend([f'$\epsilon_2 = {elem:.3f}$' for elem in eps2w_arr])
plt.title('Zavisnost prosecnog broja odbiraka m1 od greske prvog tipa')
plt.xlabel('$\epsilon_1$')
plt.ylabel('prosecan br. odbiraka')
plt.xscale('log')
plt.yscale('log')
plt.show()


plt.figure(figsize=(10,8))
for i in range(len(eps1w_arr)):
    plt.plot(eps2w_arr, m1[i, :])
plt.legend([f'$\epsilon_1 = {elem:.3f}$' for elem in eps1w_arr])
plt.title('Zavisnost prosecnog broja odbiraka m1 od greske drugog tipa')
plt.xlabel('$\epsilon_2$')
plt.ylabel('prosecan br. odbiraka')
plt.xscale('log')
#plt.yscale('log')
plt.show()

Uočavamo da grafici imaju približan oblik kao oni sa predavanja. Sada ćemo skicirati grafike i za očekivani broj odbiraka iz druge klase.

plt.figure(figsize=(10,8))
for i in range(len(eps1w_arr)):
    plt.plot(eps2w_arr, m2[i, :])
plt.legend([f'$\epsilon_1 = {elem:.3f}$' for elem in eps1w_arr])
plt.title('Zavisnost prosecnog broja odbiraka m2 od greske drugog tipa')
plt.xlabel('$\epsilon_2$')
plt.ylabel('prosecan br. odbiraka')
plt.xscale('log')
plt.yscale('log')
plt.show()

plt.figure(figsize=(10,8))
for j in range(len(eps2w_arr)):
    plt.plot(eps1w_arr, m2[:, j])
plt.legend([f'$\epsilon_2 = {elem:.3f}$' for elem in eps2w_arr])
plt.title('Zavisnost prosecnog broja odbiraka m2 od greske prvog tipa')
plt.xlabel('$\epsilon_1$')
plt.ylabel('prosecan br. odbiraka')
plt.xscale('log')
plt.yscale('log')
plt.show()

Klasifikacija

Fiksiraćemo obe greške na vrednost 0.01 i izvršiti klasifikaciju generisanih odbiraka

eps1w = 0.01
eps2w = 0.01
A = (1-eps1w) / eps2w
B = eps1w / (1-eps2w)
a = -np.log(A)
b = -np.log(B)
        
y_true = np.zeros(data.shape[1]) 
y_true[X1.shape[0]:] = 1
y_predict = np.zeros(np.shape(y_true))
start = 0 
Sm = 0

M1 = X1.mean(axis=0).reshape(-1,1)
M2 = X2.mean(axis=0).reshape(-1,1)
S1 = np.cov(X1.T)
S2 = np.cov(X2.T)

for n in range(data.shape[1]):
    f1 = gaussian_pdf(data[:,n].reshape(-1,1), M1, S1)
    f2 = gaussian_pdf(data[:,n].reshape(-1,1), M2, S2)
    h = -np.log(f1/f2)
    Sm += h

    if(Sm < a ):
        y_predict[start:n+1] = 0
        start = n+1
        Sm = 0
    elif(Sm > b):
        y_predict[start:n+1] = 1
        start = n+1
        Sm = 0

cm = confusion_matrix(y_true,y_predict)
plt.figure()
sns.heatmap(cm, annot=True, fmt='g')
plt.title('Klasifikacija pomocu Wald-ovog klasifikatora')
plt.show()

print(f'Greška prve vrste iznosi: {cm[0,1]/sum(cm[0,:]):.3f}')
print(f'Greška druge vrste iznosi: {cm[1,0]/sum(cm[1,:]):.3f}')

Greška prve vrste iznosi: 0.000
Greška druge vrste iznosi: 0.000

Zadatak 3

3.1. Generisati tri klase dvodimenzionalnih oblika. Izabrati funkciju gustine verovatnoće oblika tako da klase budu linearno separabilne.

a) Za tako generisane oblike izvšiti projektovanje linearnog klasifikatora jednom od tri iterativne procedure. Rezultate prikazati u obliku matrice konfuzije. Detaljno opisati postupak klasifikacije.

b) Ponoviti prethodni postupak korišćenjem metode željenog izlaza. Analizirati uticaj elemenata u matrici željenih izlaza na konačnu formu linearnog klasifikatora.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
import seaborn as sns
from scipy.stats import mode

Odbari parametara normalne raspodele i generisanje odbiraka

M1 = np.array([2, 2])
M2 = np.array([9, 2])
M3 = np.array([4.5, 8.5])

S = np.array([[0.8, 0], [0, 0.8]]) # Koristila sam istu kovariacionu matricu za sve 3 klase

N=100
np.random.seed(42)
X1 = np.random.multivariate_normal(M1, S, N).T
X2 = np.random.multivariate_normal(M2, S, N).T
X3 = np.random.multivariate_normal(M3, S, N).T

plt.figure(figsize=(8, 6))
plt.plot(X1[0,:], X1[1,:], 'x', label='Klasa 1')
plt.plot(X2[0,:], X2[1,:], 'x', label='Klasa 2')
plt.plot(X3[0,:], X3[1,:], 'x', label='Klasa 3')
plt.legend()
plt.title("Tri linearno separabilne klase dvodimenzionalnih oblika")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

a) Druga numerička metoda

Metoda za koju sam se odlučila je druga numerička metoda. Pošto imamo problem više kalsa, potrebno je projektovati klasifikator za svaki par klasa. U ovo slučaju ćemo ukupno imati 3 klasifikatora (1 i 2; 1 i 3; 2 i 3). Potom ćemo pomoću njih klasifikovati odbirke i to tako što ćemo odbirak smeštati u onu klasu koja je dobila najviše glasova.

Linearni klasifikator ima sledeći oblik klasifikacione funkcije: $$h(x) = V^{T}X + v_0 $$ Ukoliko je $h(x)$ veće od nule odbirak pripada drugoj, a ukoliko je manje od nule pripada prvoj klasi. Jednačine na osnovu kojih dobijamo nepoznate parametre su sledeće:$$ s = \frac{-\eta_1/\sigma_1^2}{-\eta_1/\sigma_1^2+\eta_2/\sigma_2^2} $$ $~$ $$ V = \left[s\Sigma_1+(1-s)\Sigma_2 \right]^{-1}(M_2-M_1) $$ $~$ $$ v_0 = - \frac{s\sigma_1^2V^TM_2+(1-s)\sigma_2^2V^TM_1}{s\sigma_1^2+(1-s)\sigma_2^2}$$

Pošto je nemoguće rešiti ove jednačine u zatvorenoj formi, potreno je primeniti iterativni numerički postupak. Ovaj postupak se oslanja na to da s ima opseg vrednosti od 0 do 1, kao i na to da $v_0$ može da se nađe u opsegu određenim vrednostima odbirka X koji je projektovan na vektor V, koji ćemo ovde označavati sa Y. Potrebno je zatim naći onaj par $s$ i $v_0$ koji pri klasifikaciji daje najmanju grešku. Zatim se na osnovu izabranog $s$ sračuna i vektor V.

def druga_numericka_metoda(X1, X2):
    N1 = X1.shape[1]
    N2 = X2.shape[1]
    
    M1_est = np.mean(X1, axis=1).reshape(-1,1)#reshape se radi da nam np.mean ne bi vratio kao array
    M2_est = np.mean(X2, axis=1).reshape(-1,1)
    S1_est = np.cov(X1)#zato je bilo bitno da nam dimenzija matrica bude u nacinu koji smo mi navikli, tj visoka matrica, a ne dugacka
    S2_est = np.cov(X2)    
    
    s = np.arange(0,1,10**-3)
    v0_opt_s = []#za svako s cuvamo optimalno v0 i gresku koju oni prave
    Neps_s = []
    
    for i in range(len(s)):
        V = np.linalg.inv(s[i]*S1_est + (1-s[i])*S2_est)@(M2_est-M1_est)
        # odredjivanje yj
        Y1 = V.T@X1
        Y2 = V.T@X2
        Y = np.sort(np.concatenate((Y1, Y2), axis=1)) 
        v0 = np.zeros(Y.shape[1]-1) # optimalna vrednost v0 za fiksno s
        Neps = np.zeros(Y.shape[1]-1) # Optimalno N_epsilon za fiksno s
        for j in range(Y.shape[1]-1):#PAZI, NE IDE SE DO POSLEDNJEG Y
            v0[j] = -(Y[0,j]+Y[0,j+1])/2
            Neps[j] += np.sum(Y1>-v0[j])
            Neps[j] += np.sum(Y2<-v0[j])  
        Neps_s.append(np.min(Neps))
        v0_opt_s.append(v0[np.argmin(Neps)])
    Neps_opt = min(Neps_s)
    v0_opt = v0_opt_s[np.argmin(Neps_s)]
    s_opt = s[np.argmin(Neps_s)]
    return (s_opt, v0_opt, Neps_opt, M1_est, M2_est, S1_est, S2_est)

len1 = X1.shape[1]
len2 = X2.shape[1]
len3 = X3.shape[1]

Y = np.zeros(len1+len2+len3)
Y[len1:len1+len2] = 1
Y[len1+len2:] = 2

Y1 = Y[0:len1].reshape(1,X1.shape[1])#ZASTO SU MENI OBRNUTE DIM ZA X1,X2,X3 OD NJENIH???
Y2 = Y[len1:len1+len2].reshape(1,X2.shape[1])
Y3 = Y[len1+len2:].reshape(1,X3.shape[1])

Klasifikator za klase 1 i 2

s_opt, v0_opt, Neps_opt, M1_est, M2_est, S1_est, S2_est = druga_numericka_metoda(X1, X2)
V = np.linalg.inv(s_opt*S1_est + (1-s_opt)*S2_est)@(M2_est-M1_est)
v0 = v0_opt

plt.figure()
plt.plot(X1[0,:], X1[1,:], 'rx')
plt.plot(X2[0,:], X2[1,:], 'bx')
plt.plot(X3[0,:], X3[1,:], 'gx')
plt.xlabel('obelezje1')
plt.ylabel('obelezje2')

x1 = np.arange(3,6,0.01)
x2 = -(v0+V[0]*x1)/V[1]
plt.plot(x1,x2,'m')
plt.title('Klasifikacija Klase 1 i 2')
plt.legend(['Klasa1', 'Klasa2', 'Klasa3'])
plt.show()

Kada projektujemo odbirke X na vektor V i kao granicu uzmemo $-v_0$ dobijamo slede'i grafik.

Y1 = V.T@X1
Y2 = V.T@X2
plt.figure()
plt.plot(Y1[0], np.zeros(Y1.shape[1]), 'rx')
plt.plot(Y2[0], np.zeros(Y2.shape[1]), 'bx')
plt.plot([-v0, -v0], [-1, 1])
plt.title('Klasifikacija Klase 1 i 2 - ilustracija projekcije')
plt.legend(['Klasa1', 'Klasa2'])
plt.show()

Klasifikator za klase 1 i 3

s_opt, v0_opt, Neps_opt, M1_est, M2_est, S1_est, S2_est = druga_numericka_metoda(X1, X3)
V = np.linalg.inv(s_opt*S1_est + (1-s_opt)*S2_est)@(M2_est-M1_est)
v0 = v0_opt

plt.figure()
plt.plot(X1[0,:], X1[1,:], 'rx')
plt.plot(X2[0,:], X2[1,:], 'bx')
plt.plot(X3[0,:], X3[1,:], 'gx')
plt.xlabel('obelezje1')
plt.ylabel('obelezje2')


x1 = np.arange(0,7,0.01)
x2 = -(v0+V[0]*x1)/V[1]
plt.plot(x1,x2,'m')
plt.title('Klasifikacija Klase 1 i 2')
plt.legend(['Klasa1', 'Klasa2', 'Klasa3'])
plt.show()

Y1 = V.T@X1
Y2 = V.T@X3
plt.figure()
plt.plot(Y1[0], np.zeros(Y1.shape[1]), 'rx')
plt.plot(Y2[0], np.zeros(Y2.shape[1]), 'gx')
plt.plot([-v0, -v0], [-1, 1])
plt.title('Klasifikacija Klase 1 i 3 - ilustracija projekcije')
plt.legend(['Klasa1', 'Klasa3'])
plt.show()

Klasifikator za klase 2 i 3

s_opt, v0_opt, Neps_opt, M1_est, M2_est, S1_est, S2_est = druga_numericka_metoda(X2, X3)
V = np.linalg.inv(s_opt*S1_est + (1-s_opt)*S2_est)@(M2_est-M1_est)
v0 = v0_opt

plt.figure()
plt.plot(X1[0,:], X1[1,:], 'rx')
plt.plot(X2[0,:], X2[1,:], 'bx')
plt.plot(X3[0,:], X3[1,:], 'gx')
plt.xlabel('obelezje1')
plt.ylabel('obelezje2')

x1 = np.arange(4,10,0.01)
x2 = -(v0+V[0]*x1)/V[1]
plt.plot(x1,x2,'m')
plt.title('Klasifikacija Klase 1 i 2')
plt.legend(['Klasa1', 'Klasa2', 'Klasa3'])
plt.show()

Y1 = V.T@X2
Y2 = V.T@X3
plt.figure()
plt.plot(Y1[0], np.zeros(Y1.shape[1]), 'bx')
plt.plot(Y2[0], np.zeros(Y2.shape[1]), 'gx')
plt.plot([-v0, -v0], [-1, 1])
plt.title('Klasifikacija Klase 2 i 3 - ilustracija projekcije')
plt.legend(['Klasa2', 'Klasa3'])
plt.show()

Konfuziona matrica

s12, v012, Neps_opt12, M1_est, M2_est, S1_est, S2_est = druga_numericka_metoda(X1, X2)

s13, v013, Neps_opt13, M1_est, M3_est, S1_est, S3_est = druga_numericka_metoda(X1, X3)

s23, v023, Neps_opt23, M2_est, M3_est, S2_est, S3_est = druga_numericka_metoda(X2, X3)

# Funkcija za donošenje odluke za jedan uzorak
def predict_sample(x):
    votes = {1: 0, 2: 0, 3: 0}

    V12 = np.linalg.inv(s12*S1_est + (1-s12)*S2_est) @ (M2_est - M1_est)
    if V12.T @ x > -v012:
        votes[2] += 1
    else:
        votes[1] += 1

    V13 = np.linalg.inv(s13*S1_est + (1-s13)*S3_est) @ (M3_est - M1_est)
    if V13.T @ x > -v013:
        votes[3] += 1
    else:
        votes[1] += 1

    V23 = np.linalg.inv(s23*S2_est + (1-s23)*S3_est) @ (M3_est - M2_est)
    if V23.T @ x > -v023:
        votes[3] += 1
    else:
        votes[2] += 1

    return max(votes, key=votes.get)

# Pravimo vektore stvarnih i predviđenih klasa
X_all = np.hstack([X1, X2, X3])  # sve tačke
y_true = np.array([1]*X1.shape[1] + [2]*X2.shape[1] + [3]*X3.shape[1])
y_pred = np.array([predict_sample(X_all[:, i]) for i in range(X_all.shape[1])])

cm = confusion_matrix(y_true, y_pred)
plt.figure()
sns.heatmap(cm, annot=True, fmt="d")
plt.title('Konfuziona matrica za linearni klasifikator')

Text(0.5, 1.0, 'Konfuziona matrica za linearni klasifikator')

Uočavamo 100% tačnosti, što smo i očekivali jer su klase linearno separabilne, pa možemo zaključiti da smo dobro odradili posao projekotvanja klasifikatora.

b) Metoda zeljenih izlaza

Kod ove metode se primenjuje sledeće pravilo klasifikacije: $$X \in \omega_1 \implies -V^TX-v_0>0$$ $$X \in \omega_2 \implies V^TX+v_0>0$$

Ovo se kompaknto može zapisati kao: $$W^{T}Z > 0$$ gde je $$\mathbf{W}^T = \begin{bmatrix} V_0 & V_1 & ... & V_n \\ \end{bmatrix}$$

$$X \in \omega_1 \implies \mathbf{Z} = \begin{bmatrix} -1 \\ -X \\ \end{bmatrix} $$

$$X \in \omega_1 \implies \mathbf{Z} = \begin{bmatrix} 1 \\ X \\ \end{bmatrix} $$

Nepoznate parametre dobijamo na osnovu jednačine:$$UW=\Gamma$$ koristeći pseudo inverz se dobija konačno rešenje $$W = \left(U^TU\right)^{-1}U^T\Gamma$$, gde $U$ predstavlja:$$\mathbf{U} = \begin{bmatrix} Z_1^T \\ Z_2^T \\ .\\ .\\ .\\ Z_n^T \\ \end{bmatrix}$$

$\Gamma$ predstavlja vektor željenih izlaza i tu se nalazi definisana vrednost željenog izlaza za svaki odirak. Važnijim odbircima, odnosno onim odbircima kojim želimo da pridodamo veću pažnju tokom klasifikacije dodeljujemo veću vrednost željenog izlaza $\gamma_i$

Ovde ćemo takođe kao u delu pod a) projektujemo po 1 klasifikator za svaki par klasa. Svim odbircima ćemo inicijalno pridruživati jednake željene izlaze (1), pa ukoliko uvidimo da takav klasifikator nije dovoljno dobar razmotrićemo detaljnije šta treba uraditi.

Klasifikator klase 1 i 2

x0 = np.ones((1, N*2))
x0[0,0:N] = x0[0,0:N] * -1
X = np.concatenate((-X1, X2), axis=1)
U = np.concatenate((x0, X), axis=0)

Gamma = np.ones((1, N*2)).T

W = np.linalg.inv(U@U.T)@U@Gamma
v0_12= W[0]
v1_12 = W[1]
v2_12 = W[2]

x1 = np.arange(5,6,0.01)

x2 = -(v0_12+v1_12*x1)/v2_12

plt.figure()
plt.plot(X1[0,:], X1[1,:], 'rx')
plt.plot(X2[0,:], X2[1,:], 'bx')
plt.plot(x1, x2, 'm')
plt.title("Linearni klasifikator na bazi zeljenih izlaza klase 1 i 2")
plt.legend(['Klasa1', 'Klasa2'])
plt.show()

Klasifikator kalse 1 i 3

x0 = np.ones((1, N*2)) 
x0[0,0:N] = x0[0,0:N] * -1
X = np.concatenate((-X1, X3), axis=1)
U = np.concatenate((x0, X), axis=0)

Gamma = np.ones((1, N*2)).T

W = np.linalg.inv(U@U.T)@U@Gamma
v0_13 = W[0]
v1_13 = W[1]
v2_13 = W[2]

x1 = np.arange(-2,10,0.01)

x2 = -(v0_13+v1_13*x1)/v2_13

plt.figure()
plt.plot(X1[0,:], X1[1,:], 'rx')
plt.plot(X3[0,:], X3[1,:], 'gx')
plt.plot(x1, x2, 'm')
plt.title("Linearni klasifikator na bazi zeljenih izlaza klase 1 i 3")
plt.legend(['Klasa1', 'Klasa3'])
plt.show()

Klasifikator klase 2 i 3

x0 = np.ones((1, N*2))
x0[0,0:N] = x0[0,0:N] * -1
X = np.concatenate((-X2, X3), axis=1)
U = np.concatenate((x0, X), axis=0)

Gamma = np.ones((1, N*2)).T

W = np.linalg.inv(U@U.T)@U@Gamma
v0_23 = W[0]
v1_23 = W[1]
v2_23 = W[2]

x1 = np.arange(2,10,0.01)

x2 = -(v0_23+v1_23*x1)/v2_23

plt.figure()
plt.plot(X2[0,:], X2[1,:], 'bx')
plt.plot(X3[0,:], X3[1,:], 'gx')
plt.plot(x1, x2, 'm')
plt.title("Linearni klasifikator na bazi zeljenih izlaza klase 2 i 3")
plt.legend(['Klasa2', 'Klasa3'])
plt.show()

Konfuziona matrica

N1 = X1.shape[1]
N2 = X2.shape[1]
N3 = X3.shape[1]

X = np.concatenate((X1,X2,X3), axis=1)
Y = np.zeros(X.shape[1])
Y[0:N1] = 1
Y[N1:N1+N2] = 2
Y[N1+N2:] = 3
Y_pred = np.zeros((N1+N2+N3, 3))

i1 = (X[0,:]*v1_12 + X[1,:]*v2_12 + v0_12) <= 0
i2 = (X[0,:]*v1_12 + X[1,:]*v2_12 + v0_12) > 0    
Y_pred[i1,0] = 1
Y_pred[i2,0] = 2

i1 = (X[0,:]*v1_13 + X[1,:]*v2_13 + v0_13) <= 0
i2 = (X[0,:]*v1_13 + X[1,:]*v2_13 + v0_13) > 0    
Y_pred[i1,1] = 1
Y_pred[i2,1] = 3

i1 = (X[0,:]*v1_23 + X[1,:]*v2_23 + v0_23) <= 0
i2 = (X[0,:]*v1_23 + X[1,:]*v2_23 + v0_23) > 0    
Y_pred[i1,2] = 2
Y_pred[i2,2] = 3

Y_pred_final = mode(Y_pred, axis=1).mode
Y_pred_final = Y_pred_final.reshape(-1,)

cm = confusion_matrix(Y, Y_pred_final)
plt.figure()
sns.heatmap(cm, annot=True, fmt="d")
plt.title('Konfuziona matrica za klasifikator na bazi zeljenih izlaza')
plt.show()

I ovde dobijamo klasifikaciju bez greške.

Ilustracija promene $\Gamma$

Trebalo bi da važi, da ukoliko odbircima prve klase dodamo veći željeni izlaz, klasifikator teži da njih bolje klasifikuje, pa će se klasifikaciona linija približiti drugoj klasi. Analogno važi ukoliko odbircima druge klase dodamo veću vrednos željenog izlaza.

x0 = np.ones((1, N*2))
x0[0,0:N] = x0[0,0:N] * -1
X = np.concatenate((-X1, X2), axis=1)
U = np.concatenate((x0, X), axis=0)

Gamma = np.ones((1, N*2)).T
W = np.linalg.inv(U@U.T)@U@Gamma
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1 = np.arange(2,8,0.01)
x2 = -(v0+v1*x1)/v2

Gamma2 = np.ones((1, N*2)).T
Gamma2[0:N,0] = Gamma2[0:N,0]*2 
W = np.linalg.inv(U@U.T)@U@Gamma2
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1_2 = np.arange(2,8,0.01)
x2_2 = -(v0+v1*x1)/v2

Gamma3 = np.ones((1, N*2)).T
Gamma3[0:N,0] = Gamma[0:N,0]*0.5
W = np.linalg.inv(U@U.T)@U@Gamma3
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1_3 = np.arange(2,8,0.01)
x2_3 = -(v0+v1*x1)/v2

plt.figure()
plt.plot(X1[0,:], X1[1,:], 'rx')
plt.plot(X2[0,:], X2[1,:], 'bx')
plt.plot(x1, x2, 'm')
plt.plot(x1_2, x2_2, 'm--')
plt.plot(x1_3, x2_3, 'm-.')
plt.legend(['klasa 1', 'klasa2', '$\Gamma_1$', '$\Gamma_2$', '$\Gamma_3$'])
plt.title("Linearni klasifikator na bazi zeljenih izlaza")
plt.show()

x0 = np.ones((1, N*2))
x0[0,0:N] = x0[0,0:N] * -1
X = np.concatenate((-X1, X3), axis=1)
U = np.concatenate((x0, X), axis=0)

Gamma = np.ones((1, N*2)).T
W = np.linalg.inv(U@U.T)@U@Gamma
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1 = np.arange(-2,10,0.01)
x2 = -(v0+v1*x1)/v2

Gamma2 = np.ones((1, N*2)).T
Gamma2[0:N,0] = Gamma2[0:N,0]*2 
W = np.linalg.inv(U@U.T)@U@Gamma2
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1_2 = np.arange(-2,10,0.01)
x2_2 = -(v0+v1*x1)/v2

Gamma3 = np.ones((1, N*2)).T
Gamma3[0:N,0] = Gamma[0:N,0]*0.5
W = np.linalg.inv(U@U.T)@U@Gamma3
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1_3 = np.arange(-2,10,0.01)
x2_3 = -(v0+v1*x1)/v2

plt.figure()
plt.plot(X1[0,:], X1[1,:], 'rx')
plt.plot(X3[0,:], X3[1,:], 'gx')
plt.plot(x1, x2, 'm')
plt.plot(x1_2, x2_2, 'm--')
plt.plot(x1_3, x2_3, 'm-.')
plt.legend(['klasa 1', 'klasa3', '$\Gamma_1$', '$\Gamma_2$', '$\Gamma_3$'])
plt.title("Linearni klasifikator na bazi zeljenih izlaza")
plt.show()

x0 = np.ones((1, N*2))
x0[0,0:N] = x0[0,0:N] * -1
X = np.concatenate((-X2, X3), axis=1)
U = np.concatenate((x0, X), axis=0)

Gamma = np.ones((1, N*2)).T
W = np.linalg.inv(U@U.T)@U@Gamma
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1 = np.arange(2,12,0.01)
x2 = -(v0+v1*x1)/v2

Gamma2 = np.ones((1, N*2)).T
Gamma2[0:N,0] = Gamma2[0:N,0]*2 
W = np.linalg.inv(U@U.T)@U@Gamma2
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1_2 = np.arange(2,12,0.01)
x2_2 = -(v0+v1*x1)/v2

Gamma3 = np.ones((1, N*2)).T
Gamma3[0:N,0] = Gamma[0:N,0]*0.5
W = np.linalg.inv(U@U.T)@U@Gamma3
v0 = W[0]
v1 = W[1]
v2 = W[2]
x1_3 = np.arange(2,12,0.01)
x2_3 = -(v0+v1*x1)/v2

plt.figure()
plt.plot(X2[0,:], X2[1,:], 'bx')
plt.plot(X3[0,:], X3[1,:], 'gx')
plt.plot(x1, x2, 'm')
plt.plot(x1_2, x2_2, 'm--')
plt.plot(x1_3, x2_3, 'm-.')
plt.legend(['klasa 2', 'klasa3', '$\Gamma_1$', '$\Gamma_2$', '$\Gamma_3$'])
plt.title("Linearni klasifikator na bazi zeljenih izlaza")
plt.show()

Zadatak 4

1. Generisati po $𝑁 = 500$ dvodimenzionih odbiraka iz četiri klase koje će biti linearno separabilne. Preporuka je da to budu Gausovski raspodeljeni dvodimenzioni oblici. Izabrati jednu od metoda za klasterizaciju (c mean metod, metod kvadratne dekompozicije) i primeniti je na formirane uzorke klasa. Izvršiti analizu osetljivosti izabranog algoritma na početnu klasterizaciju kao i srednji broj potrebnih iteracija. Takođe izvršiti analize slučaja kada se apriorno ne poznaje broj klasa.

2. Na odbircima iz prethodne tačke izabrati jednu od metoda klasterizacije (metod maksimalne verodostojnosti ili metod grana i granica) i primeniti je na formirane uzorke klasa. Izvršiti analizu osetljivosti izabranog algoritma na početnu klasterizaciju kao i srednji broj potrebnih iteracija. Takođe izvršiti analize slučaja kada se apriorno ne poznaje broj klasa.

3. Generisati po $N = 500$ dvodimenzionih odbiraka iz dve klase koje su nelinearno separabilne. Izabrati jednu od metoda za klasterizaciju koje su primenjive za nelinearno separabilne klase (metod kvadratne dekompozicije ili metod maksimalne verodostojnosti) i ponoviti analizu iz prethodnih tačaka.

import numpy as np
import matplotlib.pyplot as plt
import scipy

1. C mean metoda

Osobine algoritma:

• Broj klastera se mora unapred znati
• Početna klasterizacija se može izvršiti totalno stohastički
• Deli prostor bisektrisama na podprostore, ne vodi računa o kovarijansama klastera
• Ne garantuje globalni optimum
• Ne garantuje stabilnost konvergencije

Cilje je minimizovati kriterijumsku funkciju $J= tr(S_m^{-1}S_w)$. Pokazuje se da se ovaj izraz svodi na: $$J = \frac{1}{N} \sum_{i=1}^{L}\sum_{j=1}^{N_i}||X_j^{(i)} - M_i||^2$$ ukoliko se izvrši beljenje, odnosno važi $S_m = I$. Ovaj izraz će biti minimalan kada su svi odbirci smešteni u najbliži klaster po euklidskom rastojanju.

Algoritam vrši klasifikaciju na osnovu priraštaja kriterijumske funkcije:$$\Delta J(i,j,l) = ||X_i - M_j(l)||^2 - ||X_i - M_{X_i}(l)||^2$$

Koraci:

1. inicijalizacij $l=0$
2. napravimo početnu inicijalizaciju $\Omega (0)$ koja može biti potpuno nasumična
3. Procenimo vektor matematičkog očekivanja za svaku od kalsa
4. Svaki od odbiraka na osnovu priraštaja kriterijumske funckije $\Delta J(i,j,l)$ reklasifikujemo u odgovarajuću klasu - tražimo klasu koja daje najnegativniji priraštaj funkcije i tu kalsu proglašavamo novom za naš odbirak.
5. Ažuriramo l na l+1
6. Provera da li se algoritam završava tako što proverimo da li je došlo do bar jedne reklasifikacije. Ukoliko nije algoritam se završava

M1 = np.array([-5, 0]).reshape(-1,1)      
M2 = np.array([11, 1]).reshape(-1,1)
M3 = np.array([-4, 15]).reshape(-1,1)
M4 = np.array([15, 18]).reshape(-1,1)

S = np.array([[0.8, 0], [0, 0.8]])
S1 = S2 = S3 = S4 = S

N=500
np.random.seed(42)
X1 = np.random.multivariate_normal(M1.reshape(-1), S, N).T
X2 = np.random.multivariate_normal(M2.reshape(-1), S, N).T
X3 = np.random.multivariate_normal(M3.reshape(-1), S, N).T
X4 = np.random.multivariate_normal(M4.reshape(-1), S, N).T

plt.figure(figsize=(8, 6))
plt.plot(X1[0,:], X1[1,:], 'x', label='Klasa 1')
plt.plot(X2[0,:], X2[1,:], 'x', label='Klasa 2')
plt.plot(X3[0,:], X3[1,:], 'x', label='Klasa 3')
plt.plot(X4[0,:], X4[1,:], 'x', label='Klasa 4')
plt.legend()
plt.title("Cetiri linearno separabilne klase dvodimenzionih oblika")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

Pocetna klasterizacija

def pocetna_klasterizacija(klasteri, L=4):
    my_data = np.concatenate(klasteri, axis=1)
    random_class = np.random.randint(0, L, my_data.shape[1])
    
    initial_clusters = [[] for i in range(L)]
    for i, c in enumerate(random_class):
        initial_clusters[c].append(my_data[:, i].reshape(-1,1))
        
    initial_clusters = [np.concatenate(cluster, axis=1) for cluster in initial_clusters]

    return initial_clusters

klasteri = [X1, X2, X3, X4]
poc_klas = pocetna_klasterizacija(klasteri)
cnt = 1
plt.figure(figsize=(8, 6))
for klaster in poc_klas:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Pocetna klasterizacija')
plt.legend()
plt.show()

Cmean algoritam

def Cmean(klasteri, lmax=100):
    l = 1
    reklas = True  
    centroidi = [np.mean(klaster, axis=1).reshape(-1, 1) for klaster in klasteri]
    my_data = np.concatenate(klasteri, axis=1)
    klase = np.zeros(my_data.shape[1], dtype=int)
    
    while l <= lmax and reklas:
        novi_klasteri = [[] for i in range(len(centroidi))]
        reklas = False
         
        d = np.array([np.sum((my_data - centroid)**2, axis=0) for centroid in centroidi])
        nove_klase = np.argmin(d, axis=0)
        for i, c in enumerate(nove_klase):
            novi_klasteri[c].append(my_data[:, i])
            if c != klase[i]:
                reklas = True           
        klase = nove_klase.copy()
        
        klasteri = []
        pom = []
        
        for klaster in novi_klasteri:
            if klaster:
                klasteri.append(np.array(klaster).T)
            else:
                random_point = my_data[:, np.random.randint(my_data.shape[1])].reshape(-1, 1)
                pom.append(random_point)
                klasteri.append(random_point)
                
        for point in pom:
            next_point = False
            for klaster in klasteri:
                if (point in klaster) and (klaster.shape[1]>1):
                    klaster = klaster[klaster != point]
                    break
                
        
        centroidi_pom = []
        for klaster, centroid in zip(klasteri, centroidi):
             if klaster.shape[1] > 0:
                centroidi_pom.append(np.mean(klaster, axis=1).reshape(-1, 1))
             else:
                centroidi_pom.append(np.zeros_like(centroid))
        centroidi = centroidi_pom.copy()

        l += 1

    return klasteri, l

lu = 0
klasteri, l = Cmean(poc_klas, lmax=100)
cnt =1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Cmean Klasterizacija')
plt.legend()
plt.show()
lu += l

Testiranje na pocetnu klasterizaciju

Ponovićemo ovaj eksperiment još tri puta, tako ćemo videti kako Cmean radi na primeru ukupno 4 nasumične početne klasifikacije.

poc_klas = pocetna_klasterizacija(klasteri)
klasteri, l = Cmean(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Cmean Klasterizacija 2')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(klasteri)
klasteri, l = Cmean(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Cmean Klasterizacija 3')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(klasteri)
cnt = 1
klasteri, l = Cmean(poc_klas, lmax=100)
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Cmean Klasterizacija 4')
plt.legend()
plt.show()
lu += l

print(f'Srednji broj potrebnih iteracija je {lu/5}')

Srednji broj potrebnih iteracija je 3.2

U većini slučajeva Cmean daje zadovoljavajuću klasterizaciju, međutim u jednom slučaju nije data adekvatna klasterizacija, ovo ilustruje osobinu algoritma da neće uvek konvergirati ka globalno najboljem rešenju.

Bez apriornog zanja klasa.

Ovde ćemo pokazati klasifikacije kada pretpostavljamo da imamo 2, 3, 5 i 6 klasa.

Ls = [2, 3, 5, 6]

for L in Ls:
    poc_klas = pocetna_klasterizacija(klasteri, L)
    klasteri, l = Cmean(poc_klas, lmax=100)
    cnt = 1
    for klaster in klasteri:
        if klaster.shape[1] > 0:
            plt.figure(L)
            plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
            cnt += 1
    plt.title(f'Cmean Klasterizacija L={L}')
    plt.legend()
    plt.show()

2. Maximum Likelihood metoda

U ovoj metodi uvodimo pretpostavku da naši odbirci stižu iz raspodele koja se može podelovati kao zbir L gausovskih raspodela, koje su opisane aprironom verovatnoćom $P_i$, vektorom matematičkog očekivanja $M_i$ i kovarijacionom matricom $\Sigma_i$.

Potrebno je maksimizovati kriterijumsku funkciju: $$J = \max_{P_i, M_i, \Sigma_i}f(X_1, X_2, ... , X_N)$$

Metod se bazira na rešavnanju sledećih jednačina: $$q_i(X_j) = \frac{P_i \cdot f_i(X_j)}{f(X_j)}$$

$$P_i = \frac{1}{N} \sum_{j=1}^Nq_i(X_j)$$

$$M_i = \frac{1}{N_i}\sum_{j=1}^{N_i}q_i(X_j) \cdot X_j$$

$$\Sigma_i = \frac{1}{N_i} \sum_{j=1}^{N_i} q_i(X_j) \cdot (X_j-M_i)(X_j-M_i)^T$$

Mešutim ove jednačine nije moguće rešiti u zatvorenoj formi, pa se primenjuje iterativni postupak.

Koraci algoritma:

1. Inicijalizacija l=0

2. Zadamo inicijalnu klasterizaciju $\Omega(0)$ i procenimo $P_i(0)$, $M_i(0)$, $\Sigma_i(0)$

3. Na osnovu jednačine izračunamo $q_i^{(l)}(X_j)$

4. Na osnovu ostalih jednačina dobijamo vrednosti za $P_i(l)$, $M_i(l)$, $\Sigma_i(l)$

5. l = l + 1

6. Ponovo iz prve jednačine sračunamo $q_i^{(l+1)}(X_j)$

7. Tražimo $max(|q_i^{(l)}(X_j) - q_i^{(l-1)}(X_j)|)$, ukoliko je veći od $\Delta$ idemo u Korak 4, u suprotnom je kraj algoritma.

Odluka za $X_j$: $max(q_i(X_j)) = q_t(X_j) \implies X_j \in \omega_t$

def MaximumLikelihood(clusters, delta=0.001, lmax=100):
    l = 0
    max_dif = 2*delta
    my_data = np.concatenate(clusters, axis=1)

    P = [cluster.shape[1]/my_data.shape[1] for cluster in clusters]
    M = [np.mean(cluster, axis=1).reshape(-1, 1) for cluster in clusters]
    S = [np.cov(cluster) for cluster in clusters]

    fi = [[] for i in range(len(clusters))]
    for i in range(len(clusters)):
        for j in range(my_data.shape[1]):
            fi[i].append(gaussian_pdf(my_data[:,j].reshape(-1,1), M[i], S[i]))

    f = [poly_gaussian_pdf(my_data[:,i].reshape(-1,1),P,M,S) for i in range(my_data.shape[1])]

    q = [[] for i in range(len(clusters))]    
    for i in range(len(clusters)):
        for j in range(my_data.shape[1]):
            q[i].append(P[i]*fi[i][j]/f[j])

    while(max_dif > delta and l<lmax):
        P = [sum(q[i])/len(q[i]) for i in range(len(clusters))]

        M = []
        for i in range(len(clusters)):
            pom = np.zeros((2,1))
            for j in range(my_data.shape[1]):
                pom += q[i][j]*my_data[:, j].reshape(-1,1)
            pom = pom / (my_data.shape[1] * P[i])
            M.append(pom)

        S = []
        for i in range(len(clusters)):
            pom = np.zeros((2,2))
            for j in range(my_data.shape[1]):
                pom += q[i][j]*(my_data[:, j].reshape(-1,1) - M[i])@(my_data[:, j].reshape(-1,1) - M[i]).T
            pom = pom / (my_data.shape[1] * P[i])
            S.append(pom)

        l += 1

        fi = [[] for i in range(len(clusters))]
        for i in range(len(clusters)):
            for j in range(my_data.shape[1]):
                fi[i].append(gaussian_pdf(my_data[:,j].reshape(-1,1), M[i], S[i]))

        f = [poly_gaussian_pdf(my_data[:,i].reshape(-1,1),P,M,S) for i in range(my_data.shape[1])]

        q_new = [[] for i in range(len(clusters))]    
        for i in range(len(clusters)):
            for j in range(my_data.shape[1]):
                q_new[i].append(P[i]*fi[i][j]/f[j])

        substracted = []
        for i in range(len(q_new)):
            for j in range(len(q_new[i])):
                substracted.append(np.abs(q_new[i][j] - q[i][j]))

        max_dif = max(substracted)
        q = q_new.copy()
        
    result = [max(range(len(q)), key=lambda i: col[i]) for col in zip(*q)]
    final_clust = [[] for i in range(len(clusters))]
    for i, elem in enumerate(result):
        final_clust[elem].append(my_data[:, i])
    for i in range(len(final_clust)):
        final_clust[i] = np.array(final_clust[i]).T
        
    return final_clust, lmax

Analiza osetljivosti na početnu klasterizaciju

Izvršićemo 4 različite početne klasterizacije i ilustrovati finalnu klasterizaciju. Na osnovu ovoga ćemo i proceniti prosečan broj potrebnih iteracija.

lu = 0
k = [X1, X2, X3, X4]

poc_klas = pocetna_klasterizacija(k)
klasteri, l = MaximumLikelihood(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Maximum Likelihood Klasterizacija 1')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(k)
klasteri, l = MaximumLikelihood(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Maximum Likelihood Klasterizacija 2')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(k)
klasteri, l = MaximumLikelihood(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Maximum Likelihood Klasterizacija 3')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(k)
klasteri, l = MaximumLikelihood(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Maximum Likelihood Klasterizacija 4')
plt.legend()
plt.show()
lu += l

print(f'Srednji broj potrebnih iteracija je {lu/5}')

Srednji broj potrebnih iteracija je 80.0

Sada je klasterizacija uvek adekvatna, međutim potreban nam je znatno veći broj iteracija nego kod C mean algoritma

Bez apriornog znanja o broju kalsa

Ls = [2, 3]

for L in Ls:
    poc_klas = pocetna_klasterizacija(k, L)
    klasteri, l = MaximumLikelihood(poc_klas, lmax=100)
    cnt = 1
    for klaster in klasteri:
        if klaster.shape[1] > 0:
            plt.figure(L)
            plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
            cnt += 1
    plt.title(f'Maximum Likelihood Klasterizacija L={L}')
    plt.legend()
    plt.show()

3. Nelinearno separabilne klase - Kvadratna dekompozicija

Algoritam je analogan Cmean algoritmu sa jedinom razlikom što se distanca sada računa kao: $$d = \frac{1}{2} \cdot (X_i - M_j)^T \Sigma_j^{-1} (X_i - M_j) + \frac{1}{2} ln|\Sigma_j| - \frac{1}{2} ln|P_j|$$

Kod ovog algoritma potpuno stohastička inicijalizacija nije moguća, već je potrebno uneti neko predznanje u inicijalnu klasterizaciju.

def kvadratna_dekompozicija(clusters, lmax=100, output=False):

    l = 1
    flag = True  
    my_data = np.concatenate(clusters, axis=1)
    assigned_class = np.zeros(my_data.shape[1], dtype=int)
    centroids = [np.mean(cluster, axis=1).reshape(-1, 1) for cluster in clusters]
    S = [np.cov(cluster) for cluster in clusters]
    P = [cluster.shape[1]/my_data.shape[1] for cluster in clusters]
    

    
    while l <= lmax and flag:
        
        new_clusters = [[] for i in range(len(centroids))]
        flag = False
        
        d = np.zeros((len(centroids), my_data.shape[1]))
        
        for i in range(len(centroids)):
            for j in range(my_data.shape[1]):
                dist = 0.5*(my_data[:,j].reshape(-1,1)-centroids[i]).T @ np.linalg.inv(S[i]) @ (my_data[:,j].reshape(-1,1)-centroids[i]) + 0.5*np.log(np.linalg.det(S[i])) - 0.5*np.log(P[i])
                d[i,j] = dist[0,0]
                
        nxt_classes = np.argmin(d, axis=0)
        
        for i, c in enumerate(nxt_classes):
            new_clusters[c].append(my_data[:, i])
            if c != assigned_class[i]:
                flag = True           
        assigned_class = nxt_classes.copy()
        
        clusters = []
        random_points = []
        
        for elem in new_clusters:
            if elem:
                clusters.append(np.array(elem).T)
            else:
                random_point = my_data[:, np.random.randint(my_data.shape[1])].reshape(-1, 1)
                random_points.append(random_point)
                clusters.append(random_point)
                
        for point in random_points:
            next_point = False
            for cluster in clusters:
                if (point in cluster) and (cluster.shape[1]>1):
                    cluster = cluster[cluster != point]
                    break
                
        
        centroids_pom = []
        for cluster, centroid in zip(clusters, centroids):
             if cluster.shape[1] > 0:
                centroids_pom.append(np.mean(cluster, axis=1).reshape(-1, 1))
             else:
                centroids_pom.append(np.zeros_like(centroid))
        centroids = centroids_pom.copy()
        
        S_pom = []
        for cluster, sigma in zip(clusters, S):
             if cluster.shape[1] > 0:
                S_pom.append(np.cov(cluster))
             else:
                S_pom.append(np.zeros_like(sigma))
        S = S_pom.copy()
        
        if output:
            plot_clusters(clusters, title=f"Iteracija broj: {l}")
            print(f"Pretisni ENTER za narednu iteracije...")
            input()

        l += 1
        
    if output:
        print("Kraj klasterizacije")

    return clusters, l
    

N = 500

C = np.array([6,8]).reshape(2,1)
R = 4*np.random.rand(1,N)
Teta = 2 * np.pi * np.random.rand(1,N)
X11 = R*np.cos(Teta)
X12 = R*np.sin(Teta)
X1 = np.concatenate((X11,X12), axis=0) + C

C = np.array([6,8]).reshape(2,1)
R = 10 + 4 * np.random.rand(1,N)
Teta = 2*np.pi * np.random.rand(1,N)
X21 = R*np.cos(Teta)
X22 = R*np.sin(Teta)
X2 = np.concatenate((X21,X22), axis=0) + C

plt.figure()
plt.plot(X1[0, :], X1[1, :], 'x')
plt.plot(X2[0, :], X2[1, :], 'x')
plt.legend(['X1', 'X2'])
plt.title('Generisani odbirci')

k = [X1, X2]

Osetljivost na pocetnu klasterizaciju

lu = 0

poc_klas = pocetna_klasterizacija(k, L=2)
klasteri, l = kvadratna_dekompozicija(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Kvadratna dekompozicija 1')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(k, L=2)
klasteri, l = kvadratna_dekompozicija(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Kvadratna dekompozicija 2')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(k, L=2)
klasteri, l = kvadratna_dekompozicija(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Kvadratna dekompozicija 3')
plt.legend()
plt.show()
lu += l

poc_klas = pocetna_klasterizacija(k, L=2)
klasteri, l = kvadratna_dekompozicija(poc_klas, lmax=100)
cnt = 1
for klaster in klasteri:
    if klaster.shape[1] > 0:
        plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
        cnt += 1
plt.title('Kvadratna dekompozicija 4')
plt.legend()
plt.show()
lu += l

print(f'Srednji broj potrebnih iteracija je {lu/5}')

Srednji broj potrebnih iteracija je 12.0

Vidimo da kvadratna dekompozicija ne daje dobar rezultat. To je zato što je ova metoda osetljiva na početnu klasterizaciju. Izvršićemo ovaj algoritam još 10 puta da vidimo da li ćemo dobiti bolje rezultate.

Sada ćemo analizirati slučaj kada neko predznanje unesemo u inicijalnu klasterizaciju

l_sum = 0
m = 4

for i in range(m):
    original_clusters = [X1[:,100:], X2[:,100:]]
    initial_clusters = pocetna_klasterizacija(original_clusters, 2)
    initial_clusters[0] = np.concatenate((X1[:,0:100], initial_clusters[0]), axis=1)
    initial_clusters[1] = np.concatenate((X2[:,0:100], initial_clusters[1]), axis=1)
    klasteri, l = kvadratna_dekompozicija(initial_clusters)
    cnt = 0
    for klaster in klasteri:
        if klaster.shape[1] > 0:
            plt.plot(klaster[0, :], klaster[1, :], 'x', label=f"Klasa {cnt}")
            cnt += 1
    plt.title(f'Kvadratna dekompozicija {i+1}')
    plt.legend()
    plt.show()
    l_sum += l

print(f"Srednji broj potrebnih iteracija je {l_sum/m}")

Srednji broj potrebnih iteracija je 3.25

Zaključujemo da kvadratna dekompozicija dosta bolje konvergira kada se u početnu klasterizaciju uključi neko apriorno znanje.